{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LIIksEJ7fbxF"
      },
      "source": [
        "# Essential commands\n",
        "\n",
        "The following describes essential commands for interacting with the AlphaGenome\n",
        "API. It is broken into two sections: data and methods."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gcms9aHWNnqs"
      },
      "source": [
        "``` {tip}\n",
        "Open this tutorial in Google Colab for interactive viewing.\n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "executionInfo": {
          "elapsed": 359,
          "status": "ok",
          "timestamp": 1757661524606,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "iEs6z4rGe3lk"
      },
      "outputs": [],
      "source": [
        "# @title Install AlphaGenome\n",
        "\n",
        "# @markdown Run this cell to install AlphaGenome.\n",
        "from IPython.display import clear_output\n",
        "! pip install alphagenome\n",
        "clear_output()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rKyGK083Wwh7"
      },
      "source": [
        "# Imports"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "executionInfo": {
          "elapsed": 20968,
          "status": "ok",
          "timestamp": 1757661545857,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "V7MD3DBEfJwf"
      },
      "outputs": [],
      "source": [
        "from alphagenome.data import genome\n",
        "from alphagenome.models import dna_client\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "from google.colab import userdata"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dzkwq2tyfj0q"
      },
      "source": [
        "##  Data: model inputs"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3qR6e2XtW5IZ"
      },
      "source": [
        "### Genomic interval\n",
        "\n",
        "A genomic interval is specified using `genome.Interval`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "executionInfo": {
          "elapsed": 55,
          "status": "ok",
          "timestamp": 1757661546237,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "XIZHnO32W4Hn"
      },
      "outputs": [],
      "source": [
        "interval = genome.Interval(chromosome='chr1', start=1_000, end=1_010)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Qn9x1ArcXLVI"
      },
      "source": [
        "By default, these are human hg38 intervals. See the\n",
        "[FAQ](https://www.alphagenomedocs.com/faqs.html#what-are-the-reference-genome-versions-used-by-the-model)\n",
        "for more details on organisms and genome versions."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PCGNRUfHXOL1"
      },
      "source": [
        "#### Interval properties\n",
        "\n",
        "Access some handy properties of the interval:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "executionInfo": {
          "elapsed": 55,
          "status": "ok",
          "timestamp": 1757661546607,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "8bn73Lm3XL1C",
        "outputId": "2c89badb-00b8-4bb8-bb57-5b560fc5a469"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "1005"
            ]
          },
          "metadata": {},
          "execution_count": 5
        }
      ],
      "source": [
        "interval.center()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "executionInfo": {
          "elapsed": 53,
          "status": "ok",
          "timestamp": 1757661546932,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "fJVk-ocQXWhm",
        "outputId": "0d29ce5b-6586-4071-88f8-df41a7f449a3"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "10"
            ]
          },
          "metadata": {},
          "execution_count": 6
        }
      ],
      "source": [
        "interval.width"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BvlmfYgBXXig"
      },
      "source": [
        "#### Resize\n",
        "\n",
        "Use `genome.Interval.resize` to resize the interval around its center point:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "executionInfo": {
          "elapsed": 58,
          "status": "ok",
          "timestamp": 1757661547279,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "y72ZqANrXehY",
        "outputId": "31eb7150-c46c-4108-8d81-096d2feefd1e"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Interval(chromosome='chr1', start=955, end=1055, strand='.', name='')"
            ]
          },
          "metadata": {},
          "execution_count": 7
        }
      ],
      "source": [
        "interval.resize(100)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZrM4rMJDXkNF"
      },
      "source": [
        "#### Compare intervals\n",
        "\n",
        "We can also check the interval's relationship to other intervals:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "executionInfo": {
          "elapsed": 54,
          "status": "ok",
          "timestamp": 1757661547619,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "Ye04nJETXmBL"
      },
      "outputs": [],
      "source": [
        "second_interval = genome.Interval(chromosome='chr1', start=1_005, end=1_015)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "executionInfo": {
          "elapsed": 54,
          "status": "ok",
          "timestamp": 1757661548064,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "9yecEDzAXpIS",
        "outputId": "440af183-159a-4323-8900-f808c212f864"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "True"
            ]
          },
          "metadata": {},
          "execution_count": 9
        }
      ],
      "source": [
        "interval.overlaps(second_interval)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "executionInfo": {
          "elapsed": 53,
          "status": "ok",
          "timestamp": 1757661548425,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "tMN-FGXZXsqr",
        "outputId": "603d558b-43e3-45cf-e04e-736fc8080252"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "False"
            ]
          },
          "metadata": {},
          "execution_count": 10
        }
      ],
      "source": [
        "interval.contains(second_interval)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "executionInfo": {
          "elapsed": 62,
          "status": "ok",
          "timestamp": 1757661548777,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "sDjWXjJYXuPB",
        "outputId": "4780f7dc-7730-48e2-ade9-0c55bf5bf967"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Interval(chromosome='chr1', start=1005, end=1010, strand='.', name='')"
            ]
          },
          "metadata": {},
          "execution_count": 11
        }
      ],
      "source": [
        "interval.intersect(second_interval)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X0U15RKjXwZL"
      },
      "source": [
        "As a subtle point, AlphaGenome classes use 0-based indexing, meaning that the\n",
        "interval includes the base pair at the `start` position up to the base pair at\n",
        "the `end-1` position. See the\n",
        "[FAQ](https://www.alphagenomedocs.com/faqs.html#how-do-i-specify-a-genomic-region)\n",
        "for more on this topic."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dIDUCQKOX1Vj"
      },
      "source": [
        "### Genomic variant\n",
        "\n",
        "A `genome.Variant` specifies a genetic variant:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "executionInfo": {
          "elapsed": 54,
          "status": "ok",
          "timestamp": 1757661549119,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "R_D6AoKFXyBJ"
      },
      "outputs": [],
      "source": [
        "variant = genome.Variant(\n",
        "    chromosome='chr3', position=10_000, reference_bases='A', alternate_bases='C'\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "L9QcFhogX693"
      },
      "source": [
        "This variant changes the base `A` to a `C` at position 10\\_000 on chromosome 3\\.\n",
        "Note that the `position` attribute is 1-based to maintain compatibility with\n",
        "common public variant formats (see\n",
        "[FAQ](https://www.alphagenomedocs.com/faqs.html#how-do-i-define-a-variant) for\n",
        "more info.)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "l6SRhPTrYKY3"
      },
      "source": [
        "#### Insertions or deletions (indels)\n",
        "\n",
        "Variants can also be larger than a single base, such as insertions or deletions:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "executionInfo": {
          "elapsed": 56,
          "status": "ok",
          "timestamp": 1757661549453,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "PMmSYhSfX9K9"
      },
      "outputs": [],
      "source": [
        "# Insertion variant.\n",
        "variant = genome.Variant(\n",
        "    chromosome='chr3',\n",
        "    position=10_000,\n",
        "    reference_bases='T',\n",
        "    alternate_bases='CGTCAAT',\n",
        ")\n",
        "\n",
        "# Deletion variant.\n",
        "variant = genome.Variant(\n",
        "    chromosome='chr3',\n",
        "    position=10_000,\n",
        "    reference_bases='AGGGATC',\n",
        "    alternate_bases='C',\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_v9VY9kqYRoP"
      },
      "source": [
        "The sequence we pass for the `reference_bases` argument could differ from what\n",
        "is actually at that location in the hg38 reference genome. The model will insert\n",
        "whatever is passed as the reference and alternate bases into the sequence and\n",
        "make predictions on them."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "za3OTKasYYlb"
      },
      "source": [
        "#### Reference interval\n",
        "\n",
        "We can get the `genome.Interval` corresponding to the reference bases of the\n",
        "variant using `genome.Variant.reference_interval`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "executionInfo": {
          "elapsed": 62,
          "status": "ok",
          "timestamp": 1757661549825,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "UUyBaqWtYfc0",
        "outputId": "ec0c5839-7604-431d-948f-c99f629aaa64"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Interval(chromosome='chr3', start=9999, end=10000, strand='.', name='')"
            ]
          },
          "metadata": {},
          "execution_count": 14
        }
      ],
      "source": [
        "variant = genome.Variant(\n",
        "    chromosome='chr3', position=10_000, reference_bases='A', alternate_bases='T'\n",
        ")\n",
        "\n",
        "variant.reference_interval"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DATTTBK6YiEZ"
      },
      "source": [
        "A common use-case is to make predictions in a genome region around a variant,\n",
        "which involves resizing the `genome.Variant.reference_interval` to a sequence\n",
        "length compatible with AlphaGenome:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "executionInfo": {
          "elapsed": 55,
          "status": "ok",
          "timestamp": 1757661550153,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "YijroJwOYnfU",
        "outputId": "45f0e705-49ef-4f58-8c4b-2f22148cae01"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "1048576"
            ]
          },
          "metadata": {},
          "execution_count": 15
        }
      ],
      "source": [
        "input_interval = variant.reference_interval.resize(\n",
        "    dna_client.SEQUENCE_LENGTH_1MB\n",
        ")\n",
        "input_interval.width"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qGIpFPd_YmYc"
      },
      "source": [
        "#### Overlap with interval\n",
        "\n",
        "We can also check if a variant’s reference or alternate alleles overlap an\n",
        "`genome.Interval`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "executionInfo": {
          "elapsed": 54,
          "status": "ok",
          "timestamp": 1757661550470,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "HmiX9TELYxTD",
        "outputId": "34dc8c20-d66d-4a6d-c726-f05199c9b5b2"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Reference overlaps: False\n",
            "Alternative overlaps: True\n"
          ]
        }
      ],
      "source": [
        "variant = genome.Variant(\n",
        "    chromosome='chr3',\n",
        "    position=10_000,\n",
        "    reference_bases='T',\n",
        "    alternate_bases='CGTCAAT',\n",
        ")\n",
        "\n",
        "interval = genome.Interval(chromosome='chr3', start=10_005, end=10_010)\n",
        "\n",
        "print('Reference overlaps:', variant.reference_overlaps(interval))\n",
        "print('Alternative overlaps:', variant.alternate_overlaps(interval))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Kr0329VfZAZo"
      },
      "source": [
        "##  Data: model outputs\n",
        "\n",
        "### Track data\n",
        "\n",
        "\u003ca href=\"https://services.google.com/fh/files/misc/trackdata.png\"\u003e\u003cimg src=\"https://services.google.com/fh/files/misc/trackdata.png\" alt=\"anndata\" border=\"0\" height=500\u003e\u003c/a\u003e\n",
        "\n",
        "`track_data.TrackData` objects store model predictions. They have the following\n",
        "properties (using `tdata` as an example of a `track_data.TrackData` object):\n",
        "\n",
        "*   `tdata.values` store track predictions as a `numpy.ndarray` .\n",
        "*   `tdata.metadata` stores track metadata as a `pandas.DataFrame`. For each\n",
        "    track in the predicted values, there will be a corresponding row in the\n",
        "    track metadata describing its origin.\n",
        "*   `tdata.uns` contains additional unstructured metadata as a `dict`."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "I3HTL-NTZR6n"
      },
      "source": [
        "#### From scratch\n",
        "\n",
        "You can create your own `track_data.TrackData` object from scratch by specifying\n",
        "the values and metadata manually. The metadata must contain at least the columns\n",
        "name (the names of the tracks) and strand (the strands of DNA that the tracks\n",
        "are on):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "executionInfo": {
          "elapsed": 107,
          "status": "ok",
          "timestamp": 1757661550861,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "uwsbveEvZBbt"
      },
      "outputs": [],
      "source": [
        "from alphagenome.data import track_data\n",
        "\n",
        "# Array has shape (4,3) -\u003e sequence is length 4 and there are 3 tracks.\n",
        "values = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8], [9, 10, 11]]).astype(\n",
        "    np.float32\n",
        ")\n",
        "\n",
        "# We have both the positive and negative strand values for track1, while track2\n",
        "# contains unstranded data.\n",
        "metadata = pd.DataFrame({\n",
        "    'name': ['track1', 'track1', 'track2'],\n",
        "    'strand': ['+', '-', '.'],\n",
        "})\n",
        "\n",
        "tdata = track_data.TrackData(values=values, metadata=metadata)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qEc4ioZVZgR0"
      },
      "source": [
        "#### Resolution\n",
        "\n",
        "It’s also useful to specify the resolution of the tracks and the genomic\n",
        "interval that they come from, if you have this information available:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "executionInfo": {
          "elapsed": 52,
          "status": "ok",
          "timestamp": 1757661551179,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "JGnuY6g_ZaBL"
      },
      "outputs": [],
      "source": [
        "interval = genome.Interval(chromosome='chr1', start=1_000, end=1_004)\n",
        "\n",
        "tdata = track_data.TrackData(\n",
        "    values=values, metadata=metadata, resolution=1, interval=interval\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QeN8D_yGZlVB"
      },
      "source": [
        "Note that the length of the values has to match up with the interval width and\n",
        "resolution. Here is an example specifying that the values actually represent\n",
        "128bp resolution tracks (i.e., each number is a summary over 128 base pairs of\n",
        "DNA):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {
        "executionInfo": {
          "elapsed": 53,
          "status": "ok",
          "timestamp": 1757661551532,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "I-iPD_ZGZjMI"
      },
      "outputs": [],
      "source": [
        "interval = genome.Interval(chromosome='chr1', start=1_000, end=1_512)\n",
        "\n",
        "tdata = track_data.TrackData(\n",
        "    values=values, metadata=metadata, resolution=128, interval=interval\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HhVG6w41ZqEp"
      },
      "source": [
        "####  Converting between resolutions\n",
        "\n",
        "We can also interconvert between resolutions. For example, given 1bp resolution\n",
        "predictions, we can downsample the resolution (by summing adjacent values) and\n",
        "return a sequence of length 2:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {
        "executionInfo": {
          "elapsed": 54,
          "status": "ok",
          "timestamp": 1757661551868,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "Ce41uPkcZoWd",
        "outputId": "5abf98f4-2ad7-4a7c-d611-1d4bd49185c7"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[ 3.,  5.,  7.],\n",
              "       [15., 17., 19.]], dtype=float32)"
            ]
          },
          "metadata": {},
          "execution_count": 20
        }
      ],
      "source": [
        "interval = genome.Interval(chromosome='chr1', start=1_000, end=1_004)\n",
        "\n",
        "tdata = track_data.TrackData(\n",
        "    values=values, metadata=metadata, resolution=1, interval=interval\n",
        ")\n",
        "\n",
        "tdata = tdata.change_resolution(resolution=2)\n",
        "tdata.values"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "23Ehe3duZqC3"
      },
      "source": [
        "We can also upsample track data to get back to 1bp resolution and a sequence of\n",
        "length 4 by repeating values while preserving the sum:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {
        "executionInfo": {
          "elapsed": 54,
          "status": "ok",
          "timestamp": 1757661552224,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "dDi9OWWWZzq_",
        "outputId": "0745b825-26bc-4754-da3b-69cb3d0cda79"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[1.5, 2.5, 3.5],\n",
              "       [1.5, 2.5, 3.5],\n",
              "       [7.5, 8.5, 9.5],\n",
              "       [7.5, 8.5, 9.5]], dtype=float32)"
            ]
          },
          "metadata": {},
          "execution_count": 21
        }
      ],
      "source": [
        "tdata = tdata.change_resolution(resolution=1)\n",
        "tdata.values"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Rkne2GEwZ82L"
      },
      "source": [
        "####  Filtering\n",
        "\n",
        "`track_data.TrackData` objects can be filtered by the type of DNA strand the\n",
        "tracks are on:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {
        "executionInfo": {
          "elapsed": 54,
          "status": "ok",
          "timestamp": 1757661552589,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "RoIaGXLtaBbz",
        "outputId": "5e0e94bc-2216-4b56-dc22-0decfef682e5"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Positive strand tracks: ['track1']\n",
            "Negative strand tracks: ['track1']\n",
            "Unstranded tracks: ['track2']\n"
          ]
        }
      ],
      "source": [
        "print(\n",
        "    'Positive strand tracks:',\n",
        "    tdata.filter_to_positive_strand().metadata.name.values,\n",
        ")\n",
        "print(\n",
        "    'Negative strand tracks:',\n",
        "    tdata.filter_to_negative_strand().metadata.name.values,\n",
        ")\n",
        "print('Unstranded tracks:', tdata.filter_to_unstranded().metadata.name.values)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xogO-bWVaatH"
      },
      "source": [
        "#### Resizing\n",
        "\n",
        "We can resize the `track_data.TrackData` to be either smaller (by cropping):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {
        "executionInfo": {
          "elapsed": 55,
          "status": "ok",
          "timestamp": 1757661552930,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "clj536M9abBp",
        "outputId": "483d47f6-2fcd-4881-9b52-07d67bc6146d"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[3., 4., 5.],\n",
              "       [6., 7., 8.]], dtype=float32)"
            ]
          },
          "metadata": {},
          "execution_count": 23
        }
      ],
      "source": [
        "# Re-instantiating the original trackdata.\n",
        "tdata = track_data.TrackData(\n",
        "    values=values, metadata=metadata, resolution=1, interval=interval\n",
        ")\n",
        "\n",
        "# Resize from width (sequence length) of 4 down to 2.\n",
        "tdata.resize(width=2).values"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wl3aNKu-akaO"
      },
      "source": [
        "Or bigger (by padding with zeros):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "executionInfo": {
          "elapsed": 53,
          "status": "ok",
          "timestamp": 1757661553270,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "pHWZRppBaeC1",
        "outputId": "8628ab34-22fc-4572-849a-57414479ce92"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[ 0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.],\n",
              "       [ 0.,  1.,  2.],\n",
              "       [ 3.,  4.,  5.],\n",
              "       [ 6.,  7.,  8.],\n",
              "       [ 9., 10., 11.],\n",
              "       [ 0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.]], dtype=float32)"
            ]
          },
          "metadata": {},
          "execution_count": 24
        }
      ],
      "source": [
        "tdata.resize(width=8).values"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SRkggBayar2r"
      },
      "source": [
        "#### Slicing\n",
        "\n",
        "We can slice into specific positions of the `track_data.TrackData`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "executionInfo": {
          "elapsed": 78,
          "status": "ok",
          "timestamp": 1757661553613,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "yHcMdzN1amMz",
        "outputId": "29cf8f2a-826a-4a84-c2f3-930f25a2f220"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "slice by position:  [[ 6.  7.  8.]\n",
            " [ 9. 10. 11.]]\n",
            "slice by interval:  [[ 6.  7.  8.]\n",
            " [ 9. 10. 11.]]\n"
          ]
        }
      ],
      "source": [
        "# Get the final 2 positions only.\n",
        "print('slice by position: ', tdata.slice_by_positions(start=2, end=4).values)\n",
        "# Same, but using slice_interval:\n",
        "print(\n",
        "    'slice by interval: ',\n",
        "    tdata.slice_by_interval(\n",
        "        genome.Interval(chromosome='chr1', start=1_002, end=1_004)\n",
        "    ).values,\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "p1OiFpXGa-Ja"
      },
      "source": [
        "#### Subsetting tracks\n",
        "\n",
        "Subset (and reorder) to specific track names:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "executionInfo": {
          "elapsed": 53,
          "status": "ok",
          "timestamp": 1757661553946,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "_l9bKsuvaxtu",
        "outputId": "f39fc7aa-6be8-4f1f-d268-87725c53716d"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[ 0.,  1.],\n",
              "       [ 3.,  4.],\n",
              "       [ 6.,  7.],\n",
              "       [ 9., 10.]], dtype=float32)"
            ]
          },
          "metadata": {},
          "execution_count": 26
        }
      ],
      "source": [
        "# Get only tracks with the name 'track1'.\n",
        "track1_tdata = tdata.select_tracks_by_name(names='track1')\n",
        "track1_tdata.values"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bdTrjm4pbLIC"
      },
      "source": [
        "The metadata gets automatically filtered to `track1` too:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {
        "executionInfo": {
          "elapsed": 53,
          "status": "ok",
          "timestamp": 1757661554640,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "aErYz_nRbCID",
        "outputId": "476261f6-e480-4b48-9e84-21c9c7a819d2"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array(['track1', 'track1'], dtype=object)"
            ]
          },
          "metadata": {},
          "execution_count": 27
        }
      ],
      "source": [
        "track1_tdata.metadata.name.values"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dafFcwmBxXaI"
      },
      "source": [
        "#### Slicing and filtering with `[]` access.\n",
        "\n",
        "We can also apply operations using `[]` notation like numpy arrays, which generalizes the methods applied above to position and track axes."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {
        "executionInfo": {
          "elapsed": 62,
          "status": "ok",
          "timestamp": 1757661878870,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "I7J3XNY4xmUT",
        "outputId": "fcc9e421-66d4-4810-8c43-471e13c345c3"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "slice by interval and boolean track mask:  [[ 6.  7.]\n",
            " [ 9. 10.]]\n",
            "slice using python slices:  [[ 8.]\n",
            " [11.]]\n",
            "slice using slices and track names:  [[ 6.  7.]\n",
            " [ 9. 10.]]\n"
          ]
        }
      ],
      "source": [
        "interval_and_track_sliced = tdata[\n",
        "    genome.Interval(chromosome='chr1', start=1_002, end=1_004),\n",
        "    (tdata.metadata.name == 'track1').values,\n",
        "]\n",
        "print('slice by interval and boolean track mask: ', interval_and_track_sliced.values)\n",
        "\n",
        "position_and_track_sliced = tdata[2:4, 0:2]\n",
        "print('slice using python slices: ', position_and_track_sliced.values)\n",
        "\n",
        "position_and_track_sliced = tdata[2:4, ['track1']]\n",
        "print('slice using slices and track names: ', position_and_track_sliced.values)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lgprEnpIbPtR"
      },
      "source": [
        "Finally, if we pass in a stranded `genome.Interval` or leave unspecified as\n",
        "`None` when constructing a `track_data.TrackData`, we can reverse complement\n",
        "transform our track values in a strand-aware manner:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {
        "executionInfo": {
          "elapsed": 55,
          "status": "ok",
          "timestamp": 1757661554962,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "HeZXB4K3bMyJ",
        "outputId": "1db08211-b0af-4da9-af10-8737fdcc10cd"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[10.,  9., 11.],\n",
              "       [ 7.,  6.,  8.],\n",
              "       [ 4.,  3.,  5.],\n",
              "       [ 1.,  0.,  2.]], dtype=float32)"
            ]
          },
          "metadata": {},
          "execution_count": 29
        }
      ],
      "source": [
        "interval = genome.Interval(\n",
        "    chromosome='chr1', start=1_000, end=1_004, strand='+'\n",
        ")\n",
        "\n",
        "tdata = track_data.TrackData(\n",
        "    values=values, metadata=metadata, resolution=1, interval=interval\n",
        ")\n",
        "\n",
        "tdata.reverse_complement().values"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Hb8duvxibcXW"
      },
      "source": [
        "### Variant scoring output\n",
        "\n",
        "# \u003ca href=\"https://services.google.com/fh/files/misc/anndata.png\"\u003e\u003cimg src=\"https://services.google.com/fh/files/misc/anndata.png\" alt=\"anndata\" border=\"0\" height=500\u003e\u003c/a\u003e\n",
        "\n",
        "The output of variant scoring is in `anndata.AnnData` format, which is a way of\n",
        "scoring data together with annotation metadata. Originally developed in the\n",
        "single-cell RNA-seq field, `anndata.AnnData` is useful when you have metadata\n",
        "associated with an array of data.\n",
        "\n",
        "`anndata.AnnData` objects have the following properties (using `variant_scores`\n",
        "as an example `anndata.AnnData` object):\n",
        "\n",
        "*   `variant_scores.X` contains a `numpy.ndarray` containing the variant scores\n",
        "    per each gene in the region. This matrix has shape (`num_genes`,\n",
        "    `num_tracks`), where `num_tracks` is the number of output tracks in your\n",
        "    requested OutputType (such as `RNA_SEQ`, `DNASE`, etc.). Note that if you\n",
        "    did not use a gene-centric scorer, then `variant_scores.X` will have shape\n",
        "    (1, `num_tracks`), reflecting the fact that the variant has a single global\n",
        "    score and not per-gene score.\n",
        "*   `variant_scores.var` contains the track metadata as a `pandas.DataFrame`.\n",
        "    For every track in the scores (`num_genes`, `num_tracks`), there will be a\n",
        "    row in the track metadata explaining the track (its cell type, strand,\n",
        "    etc.).\n",
        "*   `variant_scores.obs` contains the gene metadata as a `pandas.DataFrame`.\n",
        "    Note that the gene metadata is None in the case of non gene-centric variant\n",
        "    scorers.\n",
        "*   `variant_scores.uns` contains some additional unstructured metadata that\n",
        "    logs the origin of the variant scores, namely:\n",
        "    *   The `genome.Variant` that was scored (variant\\_scores.uns\\[‘variant’\\])\n",
        "    *   The `genome.Interval` containing the interval\n",
        "        (variant\\_scores.uns\\[‘interval’\\])\n",
        "    *   The\n",
        "        [`variant scorer`](https://www.alphagenomedocs.com/api/models.html#variant-scorers)\n",
        "        that was used to generate the scores (variant\\_scores.uns\\[‘scorer’\\])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wxew4DE0bmLO"
      },
      "source": [
        "#### From scratch\n",
        "\n",
        "You are unlikely to need to create an `anndata.AnnData` object from scratch, but\n",
        "just for reference, here is how it would be done:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "metadata": {
        "executionInfo": {
          "elapsed": 56,
          "status": "ok",
          "timestamp": 1757661555286,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "e3jeE3HVbpof",
        "outputId": "dcaa6a1a-6716-4773-c5b2-10ea6c29162a"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/\u003cembedded stdlib\u003e/functools.py:912: ImplicitModificationWarning: Transforming to str index.\n",
            "/\u003cembedded stdlib\u003e/functools.py:912: ImplicitModificationWarning: Transforming to str index.\n"
          ]
        }
      ],
      "source": [
        "import anndata\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "\n",
        "# Creating a small matrix of variant scores (3 genes x 2 tracks).\n",
        "scores = np.array([[0.1, 0.2], [0.3, 0.4], [0.5, 0.6]])\n",
        "\n",
        "gene_metadata = pd.DataFrame({'gene_id': ['ENSG0001', 'ENSG0002', 'ENSG0003']})\n",
        "\n",
        "track_metadata = pd.DataFrame(\n",
        "    {'name': ['track1', 'track2'], 'strand': ['+', '-']}\n",
        ")\n",
        "\n",
        "variant_scores = anndata.AnnData(\n",
        "    X=scores, obs=gene_metadata, var=track_metadata\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7IHNdc_BbttW"
      },
      "source": [
        "## Methods: making predictions\n",
        "\n",
        "The main commands for making model predictions are:\n",
        "\n",
        "*   `dna_client.DnaClient.predict_sequence` to predict from a raw DNA string\n",
        "*   `dna_client.DnaClient.predict_interval` to predict from a genome interval (a\n",
        "    `genome.Interval`)\n",
        "*   `dna_client.DnaClient.predict_variant` to make predictions for ref and alt\n",
        "    sequences of a variant (a `genome.Variant` object)\n",
        "*   `dna_client.DnaClient.score_variant` to score the effects of a variant by\n",
        "    comparing ref and alt predictions.\n",
        "*   `dna_client.DnaClient.score_variants` the same as the above, but for scoring\n",
        "    a list of multiple variants."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "metadata": {
        "executionInfo": {
          "elapsed": 54,
          "status": "ok",
          "timestamp": 1757661555594,
          "user": {
            "displayName": "",
            "userId": ""
          },
          "user_tz": -60
        },
        "id": "lJ17YLKRQ6nq"
      },
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X10M3ojgbw4a"
      },
      "source": [
        "## Methods: visualization\n",
        "\n",
        "The main command for visualizing model predictions is:\n",
        "\n",
        "*   `alphagenome.visualization.plot_components.plot`, to turn a list of of\n",
        "    [plot components](https://www.alphagenomedocs.com/api/visualization.html#plot-components)\n",
        "    into a `matplotlib.figure.Figure`.\n",
        "\n",
        "See the\n",
        "[visualization basics guide](https://www.alphagenomedocs.com/visualization_library_basics.html)\n",
        "and\n",
        "[visualizing predictions tutorial](https://www.alphagenomedocs.com/colabs/visualization_modality_tour.html)\n",
        "for more details."
      ]
    }
  ],
  "metadata": {
    "colab": {
      "last_runtime": {},
      "provenance": [
        {
          "file_id": "1hJ2uMZ3sA8pu_UvSNikENLECC-X5XrR8",
          "timestamp": 1749822158925
        }
      ]
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}
